﻿using System;
using System.Windows.Forms;
using System.IO;

namespace RSProject
{
    public class funcStretch
    {
        // 2%线性拉伸
        public void twoLinearStretch(string P)
        {
            string s_data = P;
            string s_metadata = P + ".hdr";

            int[,,] tempArray = new int[10, 1000, 1000];
            byte[,,] dataArray = new byte[100, 1000, 1000];

            for (int i = 0; i < 10; i++)
                for (int j = 0; j < 1000; j++)
                    for (int k = 0; k < 1000; k++)
                        tempArray[i, j, k] = 0;

            //计算灰度直方图
            int[,] Histo = new int[7, 256];

            for (int i = 0; i < Global.head.bands; i++)
                for (int j = 0; j < Global.head.lines; j++)
                    for (int k = 0; k < Global.head.samples; k++)
                    {
                        for (int p = 0; p < 256; p++)
                        {
                            if (Global.intArray[i, j, k] == p)
                            {
                                Histo[i, p] = Histo[i, p] + 1;
                                break;
                            }
                        }
                    }

            //计算累计直方图
            int[,] TempAccuHisto = new int[7, 256];
            double[,] AccuHisto = new double[7, 256];
            for (int i = 0; i < Global.head.bands; i++)
                for (int p = 0; p < 256; p++)
                {
                    TempAccuHisto[i, p] = 0;
                    AccuHisto[i, p] = 0.00;
                }
            for (int i = 0; i < Global.head.bands; i++)
            {

                for (int p = 0; p < 256; p++)
                {
                    // TempAccuHisto[i,p]+=Histo[i,p];  为什么错了
                    for (int s = 0; s < p; s++)
                    {
                        TempAccuHisto[i, p] = TempAccuHisto[i, p] + Histo[i, s];
                    }

                    //AccuHisto[i, p] = TempAccuHisto[i, p] /360000;
                    //double xxx = TempAccuHisto[i, p] / 360000;
                    //AccuHisto[i, p] = xxx;
                }
            }

            for (int i = 0; i < Global.head.bands; i++)
                for (int p = 0; p < 256; p++)
                {
                    AccuHisto[i, p] = 1.000 * TempAccuHisto[i, p] / 360000;//真要吐血了，看了三个小时
                }

            //计算百分之二拉伸的左右区间端点值
            int[,] range = new int[7, 2];

            for (int i = 0; i < Global.head.bands; i++)
            {

                double mini = AccuHisto[i, 0] - 0.02;  //存储离0.02最近的灰度值
                if (mini < 0)
                    mini = -1 * mini;
                double maxi = AccuHisto[i, 0] - 0.98;  //存储离0.98最近的灰度值
                if (maxi < 0)
                    maxi = -1 * maxi;

                for (int p = 0; p < 256; p++)
                {
                    double a = AccuHisto[i, p] - 0.02;
                    if (a < 0)
                        a = -1 * a;
                    double b = AccuHisto[i, p] - 0.98;
                    if (b < 0)
                        b = -1 * b;
                    if (a < mini)
                    {
                        mini = a;
                        range[i, 0] = p;
                    }
                    if (b < maxi)
                    {
                        maxi = b;
                        range[i, 1] = p;
                    }
                }
            }

            //将左右端点值带入拉伸公式
            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 1; j < Global.head.lines - 1; j++)
                {
                    for (int k = 1; k < Global.head.samples - 1; k++)
                    {
                        if (Global.intArray[i, j, k] < range[i, 0] || Global.intArray[i, j, k] > range[i, 1])
                            tempArray[i, j, k] = Global.intArray[i, j, k];
                        else
                        {
                            tempArray[i, j, k] = (int)(255 * (Global.intArray[i, j, k] - range[i, 0]) / (range[i, 1] - range[i, 0]));
                        }
                    }
                }
            }

            //写入文件
            Global.data.write_metadata(s_metadata, "bsq",Global.head);
            BinaryWriter bw = new BinaryWriter(File.Open(s_data, FileMode.Create));

            for (int i = 0; i < 100; i++)
                for (int j = 0; j < 1000; j++)
                    for (int k = 0; k < 1000; k++)
                        dataArray[i, j, k] = 0;

            for (int i = 0; i < Global.head.bands; i++)
                for (int j = 0; j < Global.head.lines; j++)
                    for (int k = 0; k < Global.head.samples; k++)
                    {
                        dataArray[i, j, k] = Convert.ToByte(tempArray[i, j, k]);
                    }

            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 0; j < Global.head.lines; j++)
                {
                    for (int k = 0; k < Global.head.samples; k++)
                        bw.Write(dataArray[i, j, k]);
                }
            }
            bw.Close();
            MessageBox.Show("2%拉伸函数成功！");
        }

        // 任意拉伸
        public void LinearStretch(string P, int a, int b, int c, int d)
        {
            string s_data = P;
            string s_metadata = P + ".hdr";

            int[,,] tempArray = new int[10, 1000, 1000];
            byte[,,] dataArray = new byte[100, 1000, 1000];

            for (int i = 0; i < 10; i++)
                for (int j = 0; j < 1000; j++)
                    for (int k = 0; k < 1000; k++)
                        tempArray[i, j, k] = 0;

            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 0; j < Global.head.lines; j++)
                {
                    for (int k = 0; k < Global.head.samples; k++)
                    {
                        if (Global.intArray[i, j, k] < a || Global.intArray[i, j, k] > b)
                            tempArray[i, j, k] = Global.intArray[i, j, k];

                        else
                            tempArray[i, j, k] = (int)((d - c) * (Global.intArray[i, j, k] - a) / (b - a) + c);
                    }
                }
            }

            //写入文件
            Global.data.write_metadata(s_metadata, "bsq",Global.head);
            BinaryWriter bw = new BinaryWriter(File.Open(s_data, FileMode.Create));

            for (int i = 0; i < 100; i++)
                for (int j = 0; j < 1000; j++)
                    for (int k = 0; k < 1000; k++)
                        dataArray[i, j, k] = 0;

            for (int i = 0; i < Global.head.bands; i++)
                for (int j = 0; j < Global.head.lines; j++)
                    for (int k = 0; k < Global.head.samples; k++)
                    {
                        dataArray[i, j, k] = Convert.ToByte(tempArray[i, j, k]);
                    }

            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 0; j < Global.head.lines; j++)
                {
                    for (int k = 0; k < Global.head.samples; k++)
                        bw.Write(dataArray[i, j, k]);
                }
            }
            bw.Close();
            MessageBox.Show("拉伸函数成功！");
        }

        // 直方图均衡化      
        public void Equalization(int[,,] intdata,out int[,,] resultArray,Head head)
        {
            funcStatistics funcSt = new funcStatistics();
            int[,] frestr = new int[10, 256];
            frestr = funcSt.Frequency(intdata,head);

            int h;
            double q;
            resultArray = new int[10, 1000, 1000];
            for (int i = 0; i < head.bands; i++)
            {
                for (int j = 0; j < head.lines; j++)
                {
                    for (int k = 0; k < head.samples; k++)
                    {
                        h = 0;
                        for (int w = 1; w <= intdata[i, j, k]; w++)
                        {
                            h = h + frestr[i, w];
                        }
                        q = 255 * h / (600 * 600 * 1.00);
                        resultArray[i, j, k] = (int)q;
                    }
                }
            }
        }

        // 直方图规定化
        public void Specification(string P)
        {
            SaveFileDialog sfd = new SaveFileDialog();
            sfd.InitialDirectory = Environment.GetFolderPath(Environment.SpecialFolder.MyDocuments);
            if (sfd.ShowDialog() == DialogResult.OK)
            {
                funcRW func = new funcRW();
                funcStretch funcSt = new funcStretch();
                int[,,] intArrayTemp = new int[10, 1000, 1000];//参考图像数据
                byte[,,] dataArrayTemp = new byte[100, 1000, 1000];//参考图像二进制数据
                Head head = new Head();//参考图像头文件

                //原始图像
                int[,,] tempArray = Global.intArray;

                //读取参考图像
                Global.data.read_metadata(P + ".hdr", ref head);
                Global.data.read_data(P, ref intArrayTemp, ref dataArrayTemp, ref head);

                funcSt.Equalization(Global.intArray, out int[,,] newintarray,Global.head);
                funcSt.Equalization(intArrayTemp, out int[,,] secondintarray,head);

                int[,] h1 = new int[10, 256];
                int[,] h2 = new int[10, 256];
                int[,] first = new int[10, 256];
                int[,] second = new int[10, 256];
                int[,] final = new int[10, 256];
                int[,,] chayi = new int[10, 256, 256];
                int[,,] mark = new int[10, 1000, 1000];

                for (int i = 0; i < head.bands; i++)
                {
                    for (int m = 0; m < 256; m++)
                    {
                        for (int j = 0; j < head.lines; j++)
                        {
                            for (int k = 0; k < head.samples; k++)
                            {
                                if (newintarray[i, j, k] == m)
                                {
                                    h1[i, m]++;
                                }
                                if (secondintarray[i, j, k] == m)
                                {
                                    h2[i, m]++;
                                }
                            }
                        }
                    }//灰度直方图
                }

                int t1 = 0, t2 = 0;
                for(int h = 0; h < 10; h++)
                {
                    for (int n = 0; n < 256; n++)
                    {
                        t1 = t1 + h1[h, n];
                        first[h, n] = t1;

                        t2 = t2 + h2[h, n];
                        second[h, n] = t2;
                    }
                    t1 = 0;
                    t2 = 0;
                }//累计直方图

                for (int x = 0; x < 10; x++)
                {
                    for (int y = 0; y < 256; y++)
                    {
                        for (int z = 0; z < 256; z++)
                        {
                            chayi[x, y, z] = Math.Abs(first[x, y] - second[x, z]);
                        }
                    }
                }//两者累计直方图的差异

                int[,] jilu = new int[10, 256];
                int minnum, minvalue;
                for (int d = 0; d < 10; d++)
                {
                    for (int g = 0; g < 256; g++)
                    {
                        minnum = 0;
                        minvalue = chayi[d, g, 0];//假设1的灰度=g和2灰度=0的累计频数差异最小
                        for (int f = 0; f < 256; f++)
                        {
                            if (minvalue > chayi[d, g, f])
                            {
                                minvalue = chayi[d, g, f];//如果1的灰度=g和2的灰度=f的累计频数差异最小
                                minnum = f;//记录波段2的灰度值
                            }
                        }
                        jilu[d, g] = minnum;//g为波段一的灰度，minnum为与之差异最小的波段2的灰度
                    }
                }//记录最小差异

                for (int i = 0; i < head.bands; i++)
                {
                    for (int a4 = 0; a4 < 256; a4++)
                    {
                        for (int j = 0; j < head.lines; j++)
                        {
                            for (int k = 0; k < head.samples; k++)
                            {
                                if (mark[i, j, k] != 1)
                                {
                                    if (newintarray[i, j, k] == a4)
                                    {
                                        newintarray[i, j, k] = jilu[i, a4];
                                        mark[i, j, k] = 1;
                                    }

                                }

                            }

                        }
                    }
                }

                for (int q = 0; q < head.bands; q++)
                {
                    for (int j = 0; j < head.lines; j++)
                    {
                        for (int k = 0; k < head.samples; k++)
                        {
                            int a = newintarray[q, j, k];
                            dataArrayTemp[q, j, k] = Convert.ToByte(a);
                        }
                    }
                }

                Global.data.ToBSQ(sfd.FileName, Global.head, dataArrayTemp);
            }
        }

        // OIF
        public double CalOIF(int a, int b, int c)
        {
            ////计算7个标准差，直接从RSProject复制而来。  
            //计算相关系数矩阵，如上复制而来
            double[,] Cov = new double[Global.head.bands, Global.head.bands];
            double[] ex = new double[Global.head.bands];
            double[] ex2 = new double[Global.head.bands];
            double[] dx = new double[Global.head.bands];//这里指标准差
            double[,] Cor = new double[Global.head.bands, Global.head.bands];//相关系数

            //初始化
            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 0; j < Global.head.bands; j++)
                {
                    Cov[i, j] = 0;
                }
                ex[i] = 0;
                ex2[i] = 0;
            }

            //计算标准差
            for (int i = 0; i < Global.head.bands; i++)
            {
                for (int j = 0; j < Global.head.lines; j++)
                    for (int k = 0; k < Global.head.samples; k++)
                    {
                        ex[i] += Global.intArray[i, j, k];
                        ex2[i] += Math.Pow(Global.intArray[i, j, k], 2);
                    }
                ex[i] /= (Global.head.lines * Global.head.samples);
                ex2[i] /= (Global.head.lines * Global.head.samples);
                dx[i] = Math.Sqrt(ex2[i] - Math.Pow(ex[i], 2));
            }

            //计算协方差
            for (int m = 0; m < Global.head.bands; m++)
            {
                for (int n = 0; n < Global.head.bands; n++)
                {
                    for (int j = 0; j < Global.head.lines; j++)
                    {
                        for (int k = 0; k < Global.head.samples; k++)
                        {
                            Cov[m, n] += (Global.intArray[m, j, k] - ex[m]) * (Global.intArray[n, j, k] - ex[n]);
                        }
                    }
                    Cov[m, n] /= (Global.head.lines * Global.head.samples);
                }
            }

            //计算相关系数矩阵
            for (int i = 0; i < Global.head.bands; i++)
                for (int j = 0; j < Global.head.bands; j++)
                {
                    Cor[i, j] = Cov[i, j] / (dx[i] * dx[j]);
                }

            double midiup = dx[a + 1] + dx[b + 1] + dx[c + 1];
            if (Cor[a, b] < 0)
                Cor[a, b] = -1 * Cor[a, b];
            if (Cor[a, c] < 0)
                Cor[a, c] = -1 * Cor[a, c];
            if (Cor[c, b] < 0)
                Cor[c, b] = -1 * Cor[c, b];
            double mididown = Cor[a, b] + Cor[a, c] + Cor[b, c];
            double midi = 1.00 * midiup / mididown;

            return midi;
        }
    }
}
